function [param,glob] = setup(param,glob,options)

    %% Create state space for customer base
    minb            = glob.minb;
    maxb            = glob.maxb;

    bgrid           = nodeunif(glob.n(1),(minb)^glob.curv,(maxb)^glob.curv).^(1/glob.curv);
    bgrid0          = bgrid;

    
    %% Create idiosyncratic productivity matrix and grid
    [Pa, agrid] = setup_rouwen(param.rho_s, param.mu_s, param.sigma_s/(1 - param.rho_s^2), glob.n(2));
    Pa = Pa';
    agrid0  = agrid;

    %% Same same except for signal sapce
%     [Psig, sig_grid] = setup_rouwen(param.rho_s, 0, param.sigma_s/(1 - param.rho_s^2), glob.Nsignal);
%     Psig = Psig';
    
    %% Create finer grid for the histogram
    
    bgridf           = nodeunif(glob.nf(1),(minb)^glob.curv,(maxb)^glob.curv).^(1/glob.curv);
    Nbf              = size(bgridf, 1);

    sf               = gridmake(bgridf, agrid);
    Nsf              = size(sf, 1);
    
    glob.bgridf      = bgridf;
    glob.sf          = sf;
    glob.Nsf         = Nsf;
    
    glob.Bf = reshape(glob.sf(:,1), glob.nf(1), glob.nf(2));
    glob.Sf = reshape(glob.sf(:,2), glob.nf(1), glob.nf(2));

    
   %% Find stationary distribution across productivities
   
    a_stat = ones(size(agrid));
    a_stat = a_stat/sum(a_stat);
    a_stat = Pa'^10000*a_stat;
            
    %% Function space and nodes (fspace adds knot points for cubic splines)
    fspace          = fundef({'spli',bgrid,0,glob.spliorder(1)},...
                             {'spli',agrid,0,glob.spliorder(2)});

    sgrid           = funnode(fspace);
    s               = gridmake(sgrid);
    Ns              = size(s,1);

    %% Find actual grid given that fspace has added points for cubic spline

    %pgrid           = s(s(:,2)==s(1,2),1);
    %sgrid           = s(s(:,1)==s(1,1),2); % This is picking one of the other state and then 

    Nb              = size(bgrid,1);
    Ne              = size(agrid,1);

    glob.B = reshape(s(:,1), glob.n(1), glob.n(2));
    glob.S = reshape(s(:,2), glob.n(1), glob.n(2));
    
    %% Compute expectations matrix
    Phi             = funbas(fspace,s);
    PnK             = kron(Pa,speye(Nb));    % This is for expected value fn
    PnKf            = kron(Pa,speye(Nbf));
 
%    PnK             = kron(speye(Nb),Pa);    % This is for expected value fn
%    PnKf            = kron(speye(Nbf),Pa);


    %% Compute QZ matrix for approximating stationary distribution
    Qe              = kron(Pa,ones(Nb,1));
    
    Qef             = kron(Pa,ones(Nbf,1));

    glob.QeI             = kron(eye(size(Pa)),ones(Nbf,1));

    % construct Q matrix for iterating on productivity but not labor
    fspaceerg   = fundef({'spli',glob.bgridf,0,1});    
    Qb          = funbas(fspaceerg, sf(:,1));

    glob.QeNl     = dprod(Qef,Qb);

    
    
    %% get distribution of entrants over productivities
             
    sig_grid = linspace(min(agrid),max(agrid),glob.Nsignal);
    
    k     = 1/param.xi;
    sigma = k*min(exp(sig_grid)); % qmin is minimum gridpoint
    theta = sigma/k;
    
    P     = gpcdf(exp(sig_grid), k, sigma, theta);

    idx = find(a_stat < 1e-1);
    
    idx = idx(idx > glob.n(2)/2);
    abar = agrid(idx(1));
    
    P      =  [P(1), diff(P)];
    P = P + (1 - gpcdf(exp(sig_grid(end)), k, sigma, theta))/length(P);
    
    P(sig_grid > abar) = 0;
    P(sig_grid <= abar) = P(sig_grid <= abar) + (1 - sum(P(sig_grid <= abar)))/sum(sig_grid <= abar);
    
    glob.a_stat_E = P;
     
    %% Transition matrix from signal to realized productivity
     
    % just get a matrix that redistributes mass from the finer signal grid
    % to the coarser productivity grid
    
    fspaceerg   = fundef({'spli',agrid,0,1});
    entrant     = sig_grid;
    glob.Qshift      = funbas(fspaceerg, entrant');
    
    glob.sig_grid = sig_grid';
    
    %% evenly distribute mass across labor choices
     G               = zeros(Nsf,1);

    for i = 1:glob.n(2)
        idx = find(sf(:,2) == agrid(i));
        
        G(idx) = glob.a_stat_E(i)/length(idx);
    end

    glob.entrant_prod_dist = G;
    
    %% get distribution of entrants over both states now
        % Distribution of entrant firms
           
    %% Create other one time only basis matrices
    Phi2            = splibas(agrid0,0,glob.spliorder(2),s(:,2));
%    Phi3            = splibas(zgrid0,0,glob.spliorder(3),s(:,3));

    %% Declare additional global variables
    glob.bgrid0     = bgrid0;
    glob.bgrid      = bgrid;

    glob.sgrid0     = agrid0;
    glob.sgrid      = agrid;

%     glob.zgrid0     = zgrid0;
%     glob.zgrid      = zgrid;

    glob.Ps          = Pa;
%     glob.Pz          = Pz;

    glob.Qe         = Qe;

    glob.Ne         = Ne;
    glob.Np         = Nb;

    glob.fspace     = fspace;
    glob.s          = s;
    glob.Ns         = Ns;
    glob.Nb         = Nb;
    
    glob.Phi2       = Phi2;
%    glob.Phi3       = Phi3;

    glob.Phi        = Phi;
    glob.PnK        = PnK;
    glob.PnKf       = PnKf;
    
    glob.agrid = agrid;
    glob.a_stat = a_stat;
    
    glob.Qef        = Qef;
    
end
